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Starting from the recently proposed energy-based deviational formulation for solving the Boltzmann equation 
[J. -P. Peraud and N. G. Hadjiconstantinou, Phys. Rev. B 84, 2011], which provides significant computational 
speedup compared to standard Monte Carlo methods for small deviations from equilibrium, we show that 
additional computational benefits are possible in the limit that the governing equation can be linearized. The 
proposed method exploits the observation that under linearized conditions (small temperature differences) 
the trajectories of individual deviational particles can be decoupled and thus simulated independently; this 
leads to a particularly simple and efficient algorithm for simulating steady and transient problems in arbitrary 
three-dimensional geometries, without introducing any additional approximation. 



In a previous paper 1 , we presented a low variance 
Monte Carlo method for solving the Boltzmann trans- 
port equation (BTE) for phonons in the relaxation-time 
approximation whereby computational particles simulate 
only the deviation from an equilibrium distribution. The 
benefits of such control- variate formulations 2 , which we 
will refer to as deviational, are twofold: first, in the 
limit of small temperature differences, deviational meth- 
ods exhibit substantial computational speedup compared 
to traditional Monte Carlo methods 1 ; this speedup in- 
creases quadratically as the characteristic temperature 
difference goes to zero. Second, by simulating only the 
deviation from equilibrium, deviational methods seam- 
lessly and automatically focus the computational effort 
on regions where it is needed and can thus be used for 
solving otherwise intractable multiscale problems. In the 
present article, we show that for problems exhibiting 
sufficiently small temperature differences such that the 
BTE can be linearized, deviational computational par- 
ticles may be treated independently, thus lending them- 
selves to a simulation algorithm that is simpler, does not 
use any approximation in space or time, and, depending 
on the application of interest, can be several orders of 
magnitude faster than the one presented in Ref. 1. 

The deviational approach can be introduced by writing 
the governing equation (with no approximation) in the 
form 
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the relaxation time (lj, p and T respectively referring to 
the angular frequency, the polarization and the temper- 
ature), / = /(x, oj,p, 0, (j)) is the occupation number of 
phonon states, is the phonon-bundle group velocity 
and f%? = [exp (Huj / f kbT eq ) — 1] _1 is a Bose-Einstein dis- 
tribution at the "control" temperature T eq (/c& denotes 
Boltzmann's constant). In Ref. 1 we showed that vari- 
ance reduction is achieved by simulating only the distri- 
bution De d (D = D(w,p) is the density of states) using 



deviational particles and adding the result due to De^ 
analytically. According to (1), the scattering process is 
implemented by removing deviational particles from the 
distribution De d (i.e. the current deviational population) 
at a rate r(cj,p, T) _1 and replacing them with particles 
drawn from the distribution D(e loc — e^)/r(a;,p, T). 

For small temperature differences, the collision op- 
erator in (1) can be linearized by writing e loc — 
e^f « (T loc - T eq )de e ^ e JdT where T ioc denotes the lo- 
cal pseudotemperature 1 ' 3 . Therefore, scattered particles 
can be drawn 1 from the distribution 
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Since this distribution does not depend on (7] oc — T eq ) 
once normalized, a particle undergoing a scattering event 
can be drawn from (the normalized form of) (2) without 
knowledge of T\ oc \ energy conservation is simply ensured 
by conserving the particle. Although this formulation 
was originally introduced 1 as a means of truncating the 
discretization of a semi-infinite simulation domain (by 
limiting the region where computational cells were used), 
here, we show that this formulation can be used through- 
out the computational domain with considerable compu- 
tational benefits. By removing the need for sampling 
Ti oc before processing phonon scattering, the integra- 
tion timestep and computational cells found in standard 
Monte Carlo approaches 4 are, in fact, unnecessary. In- 
stead, the algorithm proceeds by simulating each particle 
independently and is therefore significantly simpler, re- 
quires no discretization in space and time-thereby avoid- 
ing the associated numerical error-requires significantly 
less storage and, depending on the problem of interest, 
can be several orders of magnitude more computationally 
efficient. 

The proposed algorithm for simulating a particle tra- 
jectory between t = and t = t final is as follows: 

I Draw the initial properties (sign s, position xo, fre- 
quency cjo, polarization po, direction SIq, and the 
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resulting group velocity vector V Pj o) of the parti- 
cle. For time-dependent calculations, also set up 
the initial time to of the particle (see below). 

II Calculate the traveling time until the first scat- 
tering (relaxation) event: uniformly draw a ran- 
dom number R e]0, 1[ and calculate At = 
-T(uj ,Po,T eq )\n(R) 

III Calculate x new = xo+V^oAt. Search for collisions 
with system boundaries in the time interval At. 

IVa If a collision with a system boundary occurs, say 
at Xfr, set x new = X5 and update the internal time 
tnew = to + 1 1 (x& - x ) 1 1 / 1 1 V^o 1 1 . Depending on the 
nature of the reflection (specular or diffuse), set the 
new traveling direction appropriately (as explained 
for example in Ref. 4). 

IVb If no collision with system boundaries occurs, the 
particle undergoes scattering at position x new = 
x new . The internal time is updated to t new = 
to + At. New frequency uo ne w and polarization p new 
are then drawn from (2). A new traveling direction 
is also chosen: in this work, we consider isotropic 
scattering, but this can easily be generalized to 
non-isotropic scattering. From these parameters, 
a new velocity vector V^ new can be defined. The 
particle sign remains unchanged by scattering. 

V Sample the contribution of segment [xo,x new ] to 
macroscopic properties (see below). 

VI If t new > t final, proceed to step I to begin sim- 
ulation of the next particle; otherwise, set {.}o = 
{-}newi where {.} denotes the set of all properties 
of the particle, and return to step II. 

The total number of particles processed, TV, is de- 
termined by the total amount of deviational energy in- 
volved in the phenomenon of interest divided by the 
effective energy carried by each computational parti- 
cle, E e ff. The latter is chosen such that the result- 
ing number of computational particles balances com- 
putational cost with the need for low statistical uncer- 
tainty. The contribution of initial and boundary con- 
ditions to the deviational population can be treated by 
specialized source terms. Denoting the sum of all source 
terms (including boundary and initial conditions) by 
Q(x,co>, Q,p,t), each particle's initial time to is randomly 
drawn by inverting the generalized cumulative distribu- 
tion ft, =0 J2 p J J J QdxdujdQdt' . For example, in a fi- 
nite ID system parametrized by the space coordinate 
x, the contribution of the initial condition (say initial 
temperature T^x) at t = 0) to Q is (47r)- 1 D|ef \S(t) = 
C UiP \Ti(x) -T eq \5(t), where C w , p = (Air^Dde^JdT; 
the contribution of an isothermal boundary at x = 
and at temperature T^i) to the half space x > is 
(A^D^Yg-eJ^HlYg-e,) = V g cos(0)C UiP \T b (t) - 
T eq \5(x)H[cos(9)], where is the angle with respect to 
the x > direction, and H the Heaviside step function. 



We now discuss the sampling process in more detail. 
Let I g (t') = J2 P I I J(^7r)- 1 Dge d (t f )dxdujdn be the 
macroscopic property of interest (at time t') in terms 
of a general microscopic property g = g(x,(jj,p,Q). Re- 
calling that the deviational simulation approximates the 
distribution e d in phase space using deviational (compu- 
tational) particles 1 , the estimate of I g (t') is given by 

I g (t') = £ eff J2 s i9[Mt')^i(t'),Pi(t')^i(t')} (3) 

i 

where symbols have their usual meanings and is the 
sign of deviational particle i. For example, if the quantity 
of interest is the z-component of the heat flux vector in 
some region of space R with volume /i(R) and defined by 
the characteristic function xr, then g = \ g -e z xn/ MR) 
and thus particle i only contributes to I g (t') if x^(^) [its 
position at t' — calculated by linear interpolation between 
(x ,£o) and (x new ,t new )] is in R. 

As in standard Monte Carlo methods, steady prob- 
lems can be sampled by replacing ensemble- averaging 
with time-averaging I g (ss) = (1/T) I g (t')dt' = 

(S e ff/T) ft>=t 7 ' s iddt f over a time period T, provided 
sufficient time t ss has passed for steady conditions to 
prevail. Computational benefits can be realized by not- 
ing that for steady conditions to be possible, the sys- 
tem must be under the influence of only steady particle 
sources (Q ^ Q(t)). By taking the limit T — » +00, the 
influence of initial conditions vanishes, allowing the sim- 
ulation to directly solve for — and thus focus all computa- 
tional effort on — the steady state. Particles are sampled 
over their complete trajectories, from emission (by the 
steady sources) to termination (which happens for ex- 
ample through absorption by a boundary), using 

1 9 = £ef f J2 s i I g[Mt)Mt)Mt),nMdt ( 4 ) 

i ^ 

because in this limit the effective deviational power from 
the steady sources S e ff = £ e ff/(t ss + T) reduces to 

£ejj/r. ' 

Mathematical proofs of this statement can be found 
in the linear transport theory literature (see for example 
Ref. 5). In the case of the heat flux in the z-direction 
(9 — ^zXr/M-R)) averaged over the domain R discussed 
above, equation (4) reduces to £ e ff s^//i(R) where 
Ci is the total algebraic length traveled in the z-direction 
by particle i while in R (can be negative if traveling in 
negative direction). 

The proposed algorithm has been extensively validated 
using a number of test problems 6 including the thin film 
problem described in Ref. 1 for which an analytical so- 
lution exists. Here, we present simulation results from 
two problems of practical interest (we use the same ma- 
terials and phonon properties as in Ref. 1 and 6). First, 
we consider the transient thermo-reflectance (TTR) ex- 
periment presented in Ref. 7 and used in Ref. 8 as a 
thermal conductivity spectroscopy technique. Using the 
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algorithm described above, we simulate the thermal re- 
sponse of a thin film of aluminum on a substrate of sili- 
con after a laser pulse irradiates the surface and provides 
localized heating at t = 0. More details on the prob- 
lem formulation can be found in Ref. 1, where it is also 
shown that the deviational formulation enabled the sim- 
ulation of the temperature field in this three-dimensional 
problem for several nanoseconds (due to the small tem- 
perature differences involved, simulation using standard 
Monte Carlo methods is too expensive). The additional 
speedup due to the present algorithm allows us to cal- 
culate the response to a single pulse up to 10 /is (Fig. 
1). Ultimately, we expect this improvement to be in- 
valuable towards the computational description of the 
phonon spectroscopy experiment discussed in Ref. 8 and 
9. 



■ Proposed method 
Timestep based variance-reduced 




Time (s) 

FIG. 1. Surface temperature (calculated as the spatial average 
in a cylinder of radius 10 fim and height 5nm) in the TTR 
experiment as a function of time, calculated with the variance- 
reduced Monte Carlo method using timesteps (see Ref. 1), and 
with the proposed method. The latter reaches significantly 
longer times. 

As a second application, we consider the determination 
of the thermal conductivity of complex periodic nanos- 
tructures, which has recently received a lot of attention in 
the literature 3,10 ' 11 . Here we consider a periodic nanos- 
tructure with a unit cell as shown in Fig. 2 in the pres- 
ence of a temperature gradient in the ^-direction. By 
calculating the heat flux in the direction of the gradient, 
we can determine the "effective" thermal conductivity of 
the nanostructure. Instead of considering an equilibrium 
T eq that is spatially constant, we allow the latter to vary 
in space. This approach has been shown to improve vari- 
ance reduction 12 because it allows the control tempera- 
ture to follow the physical temperature more closely; it is 
particularly convenient for imposing external fields such 
the one considered here, in which T eq (x.) varies linearly 
from Ti to T 2 > T x . With this choice of T eq , the BTE 
becomes 
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where the last term on the right hand side can be inter- 
preted as a volumetric source of deviational particles due 
to the imposed temperature gradient. When simulating 
(5), the periodic nature of the calculation is straightfor- 
wardly implemented by requiring that positive and neg- 
ative (deviational) particles individually obey periodic 
boundary conditions. Note that since the BTE is not 
linearized in (5), the source term formulation is valid for 
all deviational methods (e.g. Ref. 1). 

In order to avoid non-linearities in the response, and 
because our simulation method does not require large 
temperature differences for accuracy, we will assume 
small temperature differences (T 2 — Ti)/T <C 1, with 
To = (Ti +T2)/2; material properties, such as r(uo,p,T), 
as well as the distribution dejf jdT will be evaluated at 
To. In other words, in the linear regime, the source term 
in (5) is uniform in space. 

The simulation proceeds as outlined above (steady 
state sampling), with a few additional features due to 
the periodicity of the problem. Particles are drawn from 
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where is the polar angle (measured with respect to the 
z axis) of the particle traveling direction. Due to symme- 
try, the same number of negative and positive particles 
should be emitted. Particles exiting the domain are pe- 
riodically reinserted. 

The absence of absorbing boundaries coupled to en- 
ergy (particle) conservation results in infinitely long par- 
ticle trajectories which are impossible to track numeri- 
cally. To overcome this, we use the observation that after 
several scattering events, particle properties are almost 
completely randomized (i.e. independent from the initial 
state at emission) and thus can be terminated with only 
a small effect on the simulation accuracy. Fig. 2 illus- 
trates this for the case of the 2D nanostructure presented 
in the same figure: the mean heat flux contribution (av- 
eraged over many different particle trajectories) between 
scattering event j and j + 1, denoted by (Hj), decreases 
rapidly as j increases. Fig. 2(d) shows that after approx- 
imately 40 scattering events the statistical uncertainty in 
(Hj) (calculated in the present case using N = 8-10 6 par- 
ticles) becomes on the order of (Hj), suggesting that the 
benefit from collecting further samples is minimal and 
terminating the particle is justified. Furthermore, the 
error in the estimate of the heat flux can be controlled 
thanks to the exponential decay we observe in (Hj) af- 
ter a few scattering events. Development of a theoretical 
prediction for the number of scattering events a particle 
must undergo before it can be discarded and its depen- 
dence on the problem characteristics is the subject of on- 
going research work. For the moment, this criterion can 
be determined empirically as shown here. Figures 2(a) 
and (b) show the result obtained using this approach. 
The calculated thermal conductivities (9.8Wm _1 K _1 in 
the configuration of Fig. 2) are in agreement with previ- 
ous results 1 , while the computational time was reduced 
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by approximately 2 orders of magnitude. 
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FIG. 2. (a) Local temperature field (T — To) expressed in 
kelvin in a periodic nanostructure subject to a temperature 
gradient of — 10 6 Km -1 e z . (b) Local heat flux (Wm -2 ) in the 
z— direction, (c) Average particle contribution to the heat flux 
as a function of the particle's scattering event number. On 
average, contributions after the first scattering event amount 
to approximately 20% of the total heat flux, (d) Comparison 
between the absolute value in the heat flux contributions and 
their associated statistical uncertainty a/y/N (dashed line); 
a is the standard deviation in the heat flux as measured from 
simulation data. 

An interesting special case of the above problem is the 
calculation of thermal conductivities of thin films with 
diffuse boundaries (parallel to the z-direction) . In this 
case, using (5) and (6) to describe the imposed tempera- 
ture gradient reduces the problem dimensionality to one, 
namely the direction normal to the diffuse boundaries. 
The resulting problem is sufficiently simple (it admits 
an analytical solution) and in the case of isotropic scat- 
tering and diffuse walls is of sufficiently high symmetry, 
that the contribution to the heat flux after the first wall 
collision or the first scattering event vanishes, because 
the expected value of \V g ] • e z is zero. This observa- 
tion can be used to put in context the results presented 
in Ref. 13 where the thermal conductivity of nanostruc- 
tures was approximately calculated using a Monte Carlo 
approach which follows "test" phonon paths to their first 
free path termination (due to either a boundary or re- 
laxation). This treatment yields the correct result for 
a thin film due to the simplicity and symmetric nature 
of this problem; under more general conditions, termi- 
nating particle trajectories after the first collision event 
and assuming Fourier's law to be valid as assumed in 
Ref. 13, leads to an inaccurate answer. (In the case of 
the problem shown in Fig. 2, it leads to a value for the 
average heat flux/thermal conductivity that is approxi- 
mately 250% larger). 

Our theoretical formulation above also provides jus- 



tification for two of the assumptions used in Ref. 13, 
namely that the free paths follow a Poisson distribution 
(see equation (1)) and that "test" particles are emitted 
from any point of the nanostructure with equal proba- 
bility. The latter is only true because, as stated above, 
under linearized conditions, the source term (6) is con- 
stant (and particles can be terminated after their first 
scattering event in the thin- film problem). However, as 
stated above, for calculating the thermal conductivity of 
nanostructures, unless symmetry allows, deviational par- 
ticles need to be tracked well beyond their first free-path 
termination. This is also true for solving the Boltzmann 
equation under general conditions, since only then the 
correct non-equilibrium distribution of deviational carri- 
ers is obtained 2,12,14,15 . 

We conclude by emphasizing that the only approxima- 
tion introduced in this work comes from the assumption 
that the governing BTE can be linearized. As shown 
above, this is reasonable for a number of applications of 
interest. Under this condition, the proposed algorithm 
is in fact "more accurate" than alternative algorithms 
since it involves no timestep or spatial discretization. We 
also note that under this formulation deviational parti- 
cles share similarities with neutrons which also do not in- 
teract. Given the substantial literature on neutron trans- 
port simulation 5 , the room for improvement and gain in 
efficiency in the proposed formulation is considerable. 
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